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Abstract. 

In clades of closely related taxa, discordant genealogies due to incomplete lineage sorting 
(ILS) can complicate the detection of introgression. The D-statistic (a.k.a. the ABBA/BABA 
test) was proposed to infer introgression in the presence of ILS for a four-taxon clade. However, 
the original D-statistic cannot be directly applied to a symmetric five-taxon phylogeny, and the 
direction of introgression cannot be inferred for any tree topology. Here we explore the issues 
associated with previous methods for adapting the D-statistic to a larger tree topology, and 
propose new "Dfoil" tests to infer both the taxa involved in and the direction of introgressions 
for a symmetric five-taxon phylogeny. Using theory and simulations, we find that previous 
modifications of the D- statistic to five-taxon phylogenies incorrectly identify both the pairs of 
taxa exchanging migrants as well as the direction of introgression. The Dfoil statistics are 
shown to overcome this deficiency and to correctly determine the direction of introgressions. The 
Dfoil tests are relatively simple and computationally inexpensive to calculate, and can be easily 
applied to various phylogenomic datasets. In addition, our general approach to the problem of 
introgression detection could be adapted to larger tree topologies and other models of sequence 
evolution. 

(Keywords: introgression, incomplete lineage sorting, phylogenomic s, hybridization) 
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In phylogenomic analyses of closely related taxa, conflicting phylogenetic signals among 
loci are a common occurrence. Discordant genealogies (ones that disagree with each other and 
possibly the true species topology) represent both a challenge in determining the true species 
phylogeny and a potential source of additional information about a clade's evolutionary history 
(Maddison 1997; Degnan and Rosenberg 2009; Edwards 2009). Rapid successive speciation 
events at any time before the present can lead to discordant genealogies via incomplete lineage 
sorting (ILS), where two lineages fail to coalesce within a population, making it possible for 
either lineage to coalesce first with a less related population (Hudson 1983; Tajima 1983; Pamilo 
and Nei 1988). Discordant genealogies of closely related species can also occur through various 
forms of hybridization, ranging in scope from the introgression of alleles between existing 
species to the formation of new hybrid species (Currat et al. 2008; Twyford and Ennos 2012). 
Even though ILS and introgression/hybridization cause discordant gene trees, the difference in 
their effects on tree topologies allows them to be disentangled. While several methods for the 
detection of hybrid speciation in clades with ILS have been previously proposed (Sang and 
Zhong 2000; Meng and Kubatko 2009; Yu et al. 2012; Yu et al. 2013), here we focus on the 
detection of introgression between distinct, non-hybrid species in a clade where ILS is present. 

Several methods for distinguishing ILS from introgression in a four-taxon phylogeny 
(three species and an out-group; Fig. la) have been proposed previously. The D-statistic (a.k.a. 
the "ABBA/BABA test") looks for an imbalance in the relative frequency of two minor, 
discordant gene trees (Green et al. 2010; Durand et al. 201 1). Given a consensus phylogeny, ILS 
should produce the two minor topologies with equal frequency. However, introgression causes 
an imbalance toward a closer relationship between the two taxa exchanging alleles. Therefore, a 
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statistically significant imbalance toward one discordant topology (indicated by the allele pattern 
ABBA or BABA) indicates that introgression has occurred. 

The 4sp algorithm (Garrigan et al. 2012) also uses the relative frequencies of allelic state 
patterns to determine regions of introgression. This algorithm uses a maximum-likelihood 
approach to estimate global parameters of the species tree, then calculates a local likelihood of 
introgression for each region of the genome (Garrigan et al. 2012). Joly et al. (2009; see also 
Joly 2012) use a similar method to look for loci that have sequences with unusually low genetic 
distances. A significantly lower minimum sequence distance for a pair of taxa indicates 
introgression at specific loci. Both of these methods, and the D-statistic, allow for fine-scale 
analysis of localized introgression when genomic spatial information is available (i.e., large 
chromosomal alignments). 

When expanding phylogenetic tests for ILS and introgression beyond the simple four- 
taxon case, there are several challenges. In the presence of ILS, a four-taxon, rooted phylogeny 
(Fig. la) has only three different possible gene tree topologies. For a five-taxon species 
phylogeny (Fig. lb) in the presence of ILS there are fifteen possible gene tree topologies both 
symmetric (Fig. lc) and asymmetric (Fig. Id). In addition to more possible gene trees, there are 
an increased number of possible introgression donor-recipient pairs, and the probability 
distribution of possible gene trees must also be considered (Degnan and Salter 2005; Degnan and 
Rosenberg 2006; Twyford and Ennos 2012). 

Here, we propose a new set of statistical measures (the Dfoil tests) that not only 
overcome the inherent challenges of additional phylogenetic information in five-taxon symmetric 
phytogenies, but also leverage this added data to infer both the taxa involved in and direction of 
introgression. We also show how these statistics can be used to infer the occurrence of 
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introgression along ancestral branches, and present results that suggest a previous modification 
of the D-statistic (the "Partitioned" D-statistic; Eaton and Ree 2013) imprecisely identifies the 
taxa involved in introgression. The dynamics of all these statistics are shown theoretically and 
by simulation, followed by a discussion of their applications to data. 

Materials and Methods 

The Four-taxon D-statistic 

In order to approach testing introgression in a five-taxon phylogeny, we first describe the 
four-taxon case. (Note that in much of the literature on the multispecies coalescent these are 
referred to as three-taxon and four-taxon trees, where the outgroup is not counted. For 
consistency with previous work on the D-statistic, we do not use this terminology here.) The 
four-taxon D-statistic for introgression was formalized to test ancestral admixture between 
human and Neanderthal populations (Green et al. 2010; Durand et al. 2011). This statistic 
applies to a four-taxon asymmetric phylogeny with three in-group taxa and an out-group, 
denoted {{{P\,Pi),P?,),0) (Fig. la). All sites considered in the alignment of sequences from these 
taxa must be either mono- or biallelic, with the out-group defining the ancestral state (always 
named A) relative to the derived state (named B). Allelic state patterns for a given position in the 
alignment are given in the order P1P2P3O (e.g., ABBA; Fig. 1 and 2). Site pattern counts (e.g., 
Rabba) are the raw counts of these site types in a given region of the sequence alignment. The 
model generally assumes 0 or 1 substitutions at each site over the whole phylogeny, with a 
negligible number of reverse and convergent substitutions. This model also assumes 0 or 1 
introgressions in a region. The true is tree is supported by the patterns BBAA and BBBA, and 
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polyphyletic appearance of B in the discordant site patterns ABBA and BABA is attributed to 
either ILS or introgression (or both). 

The D-statistic is calculated as (Green et al. 2010; Durand et al. 201 1): 

_ ^ABBA — ^BABA 
^ABBAA + ^BABA 

Under ILS and no introgression, the discordant tree patterns ABBA and BABA should occur with 
equal frequency (Hudson 1983; Tajima 1983; Pamilo and Nei 1988). Introgression between Pi 
and either Pi or Pi will disproportionately increase the frequency of BABA or ABBA, 
respectively, since the introgressed pair of taxa should have relatively more shared derived (B) 
states. Therefore, the four-taxon D-statistic is a measure of inequality in the prevalence of site 
patterns that support the two possible discordant gene tree topologies. Positive values of D 
indicate Pi^P?, and negative values indicate P\^P?,, with values not differing significantly from 
zero not supporting either introgression (where <=> denotes introgression of indeterminate 
direction and => denotes a polarized introgression). This approach is not able to detect 
introgression between the two sister taxa, Pi and Pi. 



Considerations for Extending the D- statistic Beyond Four Taxa 

As the number of taxa in a phylogeny increases, introgression testing becomes 
increasingly complex due to three factors: (1) the number of introgression donor-recipient pairs 
increases geometrically, (2) the number of possible gene tree topologies also increases 
geometrically, and (3) the probability distributions of discordant gene tree topologies under ILS 
become more complex with both the size and shape of the phylogeny (Rosenberg 2002; Degnan 
and Salter 2005; Degnan and Rosenberg 2006; Rosenberg 2007). Another practical concern is 
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that the four-taxon D-statistic can only provide two answers (positive and or negative), 
corresponding to the introgressions Pi<=>P3 and P2<?$Pi. A single D- statistic is sufficient for the 
four-taxon case, but a larger phylogeny will require a system of multiple tests to distinguish 
among the (rapidly) increasing number of possible introgressions. 

There are multiple perspectives from which the D-statistic can be understood. In all 
interpretations, the D-statistic has the null hypothesis that D=0 with or without ILS, and with no 
introgression. From the perspective of site patterns, the D-statistic states that the biallelic site 
patterns ABBA and BABA should be sampled with equal frequency under the null hypothesis (Eq. 
1). The site patterns ABBA and BABA are a proxy for the two discordant gene trees 
(((P2,P3),Pi),0) and (((P\,Pi),P2),0), respectively (Fig. 2a, lines ii and iii), each of which is 
expected to occur with equal frequency. More generally, we can say that the D-statistic 
compares two sets of gene trees — the "left" and "right" terms of the numerator — whose sampling 
probabilities are expected to be equal given the probability distribution of gene trees for a 
particular species phylogeny. In a four-taxon phylogeny, the two sets of gene trees are each 
represented by only a single discordant gene tree (of the two possible discordant topologies). 
Since the two discordant gene trees are equally likely to be sampled under ILS, there is an equal 
probability of inferring that P3 is more closely related to Pi or Pi, and D equals 0 under the null 
hypothesis of no introgression. Importantly, in the D-statistic it is not necessary to calculate the 
exact probability of sampling either discordant gene tree. Because both discordant trees are due 
to ILS on the same ancestral branch, both are expected to have equal relative probability 
regardless of what the absolute probabilities are (Fig. 2a). This feature will be important in 
designing D statistics for five-taxon phylogenies. 
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Another important aspect of the four-taxon D-statistic is that it only uses shared derived 
states to infer introgression. However, introgression is equally capable of transferring the 
ancestral state (A) or the derived state (B). For example, consider a four-taxon phylogeny where 
a substitution A— >B has occurred on the internal branch P12 (Fig. la). This would ordinarily lead 
to the site pattern BBAA. However, if P3=> P\ introgression occurs, then the A state is transferred 
resulting in pattern BAAA. This means that when introgression is considered, both ABBA and its 

"inverse pattern," BAAA, offer evidence of P30P1 introgression, and the same is true for 

BABA/ABAA and P30P2. As long as both terms in the numerator of D use these inverse site 

counts, the null hypothesis of D=0 is maintained. However, the inclusion of patterns with a 
single derived state (i.e. one B) may cause complications in some types of sequence data and in 
some biological scenarios where substitution rates are unequal among terminal branches of the 
phylogeny (see Discussion). 

From these considerations of the four-taxon D-statistic, we can derive four general 
principles that can be used to test introgression in a five-taxon phylogeny. First, a system of 
multiple D-statistics will be required to distinguish among the larger number of possible donor- 
recipient combinations of introgression. Second, rather than designing a solution particular to 
this tree topology as a whole, we can discretize introgression testing into a system of taxon-by- 
taxon D-statistics by examining the relative relationships of a given taxon against two other 
(appropriately selected) taxa. Third, to maintain the null expectation of D=0, these two relative 
phylogenetic relationships must have an equal sampling probability across the distribution of all 
possible gene trees. As with the four-taxon statistic, it will not be necessary to calculate the 
actual probability values as long as gene trees can be selected in equally probable pairs (as 
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illustrated in Figure 2). Finally, inverse patterns (e.g. BABA/ABAA) both indicate the same 
potential introgressions. 
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The Five-Taxon "Partitioned" D- statistics 

We also use these principles to examine the "Partitioned" D-statistics previously 
proposed to infer inter-group introgression (Fig. 3a) in a symmetric five-taxon phylogeny (Eaton 
and Ree 2013). These statistics are: 

(2) 
(3) 

'12 ~ ~ — (4) 

Positive values of Di indicate P2<^P?, introgression, while negative values indicate Pi<^>Pi 
(Eaton and Ree 2013). Values not differing significantly from 0 indicate no introgression 
involving P3. For D2, positive and negative values indicate Pi<^Pa and P\<^Pa, respectively, 
while D2 = 0 indicates no introgressions involving Pa. In addition, the Dn statistic was proposed 
as a means to determine the introgression donor and recipient taxa. Dn assumes that if both P3 
and P\ exhibit the derived state, then the substitution must have occurred on the ancestral P34 
branch (Fig. lb). So Pi or Pa is assumed to be the introgression donor, and significant positive 
or negative values of Dn indicate the recipient is Pi or Pz, respectively. Values not deviating 
significantly from 0 indicate Pi or Pi as the introgression donor and P3 or Pa as the recipient. 



Simulations 
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We tested the new Dfoil statistics described below on a set of simulated chromosomes 
divided into 10 kb regions, each of which was allowed to evolve over a common phylogeny (Fig. 
3b and 3c). All simulated chromosomes were generated using the Python libraries for the 
EggLib coalescent simulation engine (v2.1.7; De Mita and Siol 2012). A simulated population 
size of /V e =10 6 and per-base-per-generation mutation rate of //=10" 9 were used. Introgression was 
simulated by a temporary period of migration from the introgression donor population to the 
recipient population at a fixed rate (500 individuals per generation). Biallelic site patterns were 
counted for each region. To simulate the possible effects of convergent mutations in a finite sites 
mode, five additional substitutions were added randomly to each 10 kb region at sites with one or 
more existing derived states. The four Dfoil and three Partitioned D-statistics were then 
calculated from the site pattern counts for each 10 kb window. Custom Python scripts were used 
for tabulation of site patterns and calculation of Partitioned D-statistics and Dfoil. Plots were 
generated with the matplotlib library (http://www.matplotlib.org). The simulation generation 
script and site count data sets are available through Dryad (URL PENDING). The DFOIL 
program used for all calculations and plots is for available for public use 
(http://www.bitbucket.org/jbpease/dfoil). 

Results 

The Dfoil Test for a Symmetric Five-taxon Phylogeny 

Our model describes a clade of five taxa connected by a symmetric phylogeny, denoted 
(({P\,Pi),{P?>,P*)),0), with the in-group taxa arranged in two sub-pairs (P1/P2 and P3/P4) and an 
out-group taxon (O; Fig. lb). We define that P3 and Pa diverged (at time-before-present Ti) 
before Pi and Pi (at Ti), and also that the two sub-pair lineages diverged at 7i. The labeling of 
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the taxa (Pi— Pa) is arbitrary, as long as the sub-pairings are correct and the relationships between 
the three speciation times-before-present adhere to the relationships Ti > Ti > T3 > 0. 

As in the four-taxon case, we only sample mono- or biallelic sites (with the out-group 
always represented as A). 0 or 1 introgressions are allowed per region (at time-before -present r). 
Reverse and convergent substitutions are expected to occur in negligible amounts, and are 
expected to affect all topologies equally (Durand et al. 201 1). We refer to introgressions 
occurring between one taxon from each sub-pair as an inter-group introgression (of which there 
are n = 8 possible pairings; Fig. 3a) and those between taxa in the same sub-pair as an intra- 
group introgression in = 4; Fig. 3a). Additionally, introgression between the ancestral branch Pn 
and P3 or P\ is possible when Ti>x> T3, which will be referred to as an ancestral introgression 
(n = 4; Fig. 3a). 

We propose a system of four Z)-statistics, collectively named Dfoil, to distinguish among 
the 16 possible introgressions in a symmetric five-taxon phylogeny. The name Dfoil borrows 
from the "FOIL method," a grade school mnemonic for multiplying two binomials ("First, Outer, 
Inner, Last"). We apply these labels to the four in-group taxa, and name the four Dfoil statistics 
Dfo ("first'WWs vs. "outer" =Pi/P 4 ), Dil ("inner"=/VP 3 vs. "last"=P 2 // 5 4), Dpi ("first" vs. 
"inner"), and Dol ("outer" vs. "last"). These are defined as: 



O^BABAaF ^BBBAA F ^ ABABA F ^AAABa) O^BAABA F ^BBABA F W ABBAA F ^AABAa) 

(^babaaF ^bbbaa "F ^ababa F ^aaaba) "F (^baaba F ^bbaba F ^abbaa F JIaabaa) 



„ O^ABBAaF ^BBBAA F 71-BAABA + ^AAABa) O^ABABA F J^BBABA F ^BABAA F ^aABAa) ,,r\ 

IL = 7 S — 7 ^ ' ^ 

C^ABBAaF ^BBBAA F ^-BAABA F ^aAABaJ F (tl ABABA F Tl BBABA + Tl-BABAA F W MBM J 



(^BABAaF ftfiABBA F JT-ABABA F ^aBAAa) O^ABBAA F ^ABBEA F J^BAABA F JI-bAAAa) 

(^babaaF ^babba F JT-ababa F ^abaaa) F (^abbaa F ^abbba F 71-baaba F TL baaaa ) 
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^ _ (^BAABA"!" ^BABBA + ^ABBAA ^ABAAa) (^ABABA "I" ^ABBBA ^BABAA + ^BAAAa) ^ 
(^BAABA"!" ^BABBA + ^ABBAA + ^ABAAa) "I" (^ABABA + ^ABBBA + ^BABAA + ^BAAAa) 



Each of these equations calculates the D-statistic for one of the four in-group taxa using the 
principles described previously. For example, Dfo tests Pi and describes the relative support for 
Pi being more closely related to P3 or Pa (i.e., the two taxa from the opposite sub-pair; Fig lb). 
These two relationships are inferred by sampling two sets of gene trees (indicated by biallelic 
site patterns) that have an equal total probability of being sampled under the null hypothesis 
given the distribution of gene trees for a symmetric five-taxon phylogeny ("+Dfo" and "—Dfo" in 
Fig. 2b). Under ILS alone, the relative apparent strength of the two relationships represented on 
either side of the numerator of Dfo should be equal. Therefore, Pi<=>P3 introgression will lead to 
more sampling of sites that support a closer relationship between Pi and P3, and therefore shift 
towards Dfo> 0. Alternatively, Pi<=>P4 introgression will shift towards Dfo< 0. We also apply 
the principle of inverse patterns (e.g. BABBA and ABAAA), and so both terms of all Dfoil tests 
include two pairs of inverse patterns. Dil, Dfi, and Dol are calculated identically to Dfo, with 
Pi, P3, and Pa, respectively, as the focal taxon instead of Pi, and P1/P2 used for comparison in 
Dfi and Dol. For all four tests, significant positive or negative values support a hypothesis of 
introgression for the focal taxon with one of the two taxa from the other sub-pair. 

Importantly, it should be noted that the constraint of equal probability for each of the two 
terms in the numerators of the Dfoil statistics means that an analogous set of tests for an 
asymmetric five-taxon phylogeny cannot be constructed. In the asymmetric phylogeny 
((((Pi,P2),P3),Pa)0) the relationship between Pi and P3 is not expected to be equal to Pi and Pa 
under the null. This occurs because a closer relationship between Pi and P3 only requires ILS to 
have occurred on one branch, while a closer relationship between Pi and Pa requires ILS on two 
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ancestral branches. This highlights the special property that the four-taxon asymmetric and five- 
taxon symmetric phylogenies share with respect to ILS. In both topologies, all discordant 
coalescences must occur in the root, and therefore the exact probability distributions of 
discordant topologies are not needed to construct a test of the null hypothesis of no introgression. 

Direction of Introgression and Significance 

As described thus far, the four Dfoil statistics are simply individual applications of the D- 
statistic to each of the four in-group taxa. However, all four Dfoil statistics considered 
collectively contain more information than the sum of the individual D-tests. In addition to 
identifying the taxa involved in introgression, the Dfoil statistics can also provide information 
about the direction of inter-group introgressions, specifically identifying both the donor and 
recipient taxa. When an inter-group introgression occurs in a symmetric five-taxon phylogeny, 
the relationship changes not only between the donor and recipient taxa change, but also the 
recipient taxon and the donor's sister taxon. For example, if P?,=>P\ occurs in an otherwise 
concordant gene tree, then the resultant topology becomes {{{P\,P^),Pa),Pi),0) (Fig 2b, line viii). 
The shared history of P3 and Pa means that Pa also changes its relationship to Pi. Conversely, if 
Pi=>P3 occurred the resultant topology is (((Pi,P?>),P2),Pa),0) (Fig 2b, line vii), and now the 
relationship of Pi with P3 and Pa is altered by association. Even though both relationships 
change, the introgressing taxon should change more strongly than its sister taxon, and this 
disparity informs the overall assignment of donor and recipient. In this way, the signs (+, -, or 
0) of the four Dfoil tests collectively form a signature that provides information beyond the sum 
of its individual components. The Dfoil tests individually indicate which taxa are introgressing, 
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but collectively can also identify the donor and recipient taxa (or ancestral lineage) for a given 
introgression. 

Each Dfoil tests is separately assessed to be significantly positive, significantly negative, 
or not different than 0 by a two-tailed binomial test of the raw site counts for the left and right 
terms with a 50:50 expectation. Statistically significant positive and negative values are assigned 
a sign of "+" or while non-significant values are "0." The four tests in the order {Dfo, Dil, 
Dfi, Dol) form the Dfoil signature (Table 1). The appropriate P- value cutoffs for significance 
are determined by simulation, using the tested topology with no introgression. The 0.05 tail of 
simulated components of Dfoil from regions with no introgression are used as the empirical P- 
value cutoffs. Individual P-values for each Dfoil component can be calculated; however, it is 
expected that Dfo and Dil will have approximately equal significance thresholds, as will Dfi and 
Dol. The difference in the thresholds between these two pairs of statistics should be 
approximately proportional to the difference in divergence times Ti and T3. 

Application and accuracy of the Dfoil method 

To demonstrate the effectiveness of the Dfoil method, we simulated the evolution of 
chromosomes over a common phylogeny (Fig. 3a and 3b). Along the chromosome, we 
simulated sixteen sets of n=500 loci of length 10 kb for each of the sixteen possible 
introgressions (eight inter- group, four ancestral branch, four intra- group). Between each 
introgressed region, n=100 regions of no-introgression control were included for visual contrast 
(Fig. 3d, shaded). In all sixteen cases and the control, the Dfoil signature (signs of the four Dfoil 
tests) that was observed in the simulations matched the theoretical expectations (Fig. 3d and 
Table 1). For example, the introgression Pi=>Pi (Fig. 3d, leftmost region) shows positive values 
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for Dfo, Dil, and Dfi, while Dol ~ 0, in agreement with the expected +/ + / + / 0 signature. In 
each of the eight inter-group introgressions (Fig. 3d, leftmost eight regions), a unique Dfoil 
signature clearly distinguishes each case. Ancestral introgressions involving branch Pn (Fig. 3c 
and middle of 3d) were distinguishable from cases of introgression with P3 or Pa, but the 
direction of these introgression cannot be determined by the Dfoil signature alone (see 
Discussion). As expected, the average value for all four Dfoil tests in intra-group introgressions 
(Fig. 3d, rightmost four regions) is 0. Therefore, the Dfoil tests can distinguish the taxa involved 
in inter-group and ancestral introgressions, and additionally the direction of introgression for 
inter- group introgressions. 

The Five-Taxon "Partitioned" D- statistics 

We also used the D-statistic framework described above to examine the Partitioned D- 
statistics (Eq. 2-4; Eaton and Ree 2013). As previously noted, inverse pairs of biallelic patterns 
(e.g. ABBAA/BAABA) both indicate the same underlying gene tree when introgression is 
considered. This means that the site pattern counts used in D\ and Di of the Partitioned D- 
statistic are not necessarily unique to the introgressions they propose to test. Specifically, the left 
term of Di (n ABBAA ) and right term of Di (h B aaba) are directly related, as are the right term of D\ 
(«babaa) and left term of Di (h AB aba)- Therefore, an inter-group introgression between any two 
taxa should change both Di and D2 in opposite directions. 

This inverse relationship between Di and Di means that introgression between one pair of 
taxa creates a "mirror effect" that makes the other pair of taxa falsely exhibit evidence of 
introgression. Therefore, while D\ and Di in this form can detect the occurrence of an inter- 
group introgression, the mirror effect makes it unclear which two taxa are actually exchanging 
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alleles. This effect is a consequence of using biallelic patterns where two (out of four) taxa 
exclusively share one state, since, by default, the other two taxa share the opposite state. 

In all cases of inter-group introgression, the simulation data showed an inverse 
relationship between Di and D2, confirming the mirror effect (Fig. 4). For example, the 
introgressions Pi=>Ps and P?,=$P\ (Fig. 4, left side) are expected to be Di < 0 and D2 ~ 0. 
However, in both cases Di showed positive values that mirrored the negative values of Di, 
making it falsely seem as though P2<^Pa was also occurring in both cases. A similar result was 
observed for the introgressions Pi=>Pa and Pa=>Pi (Fig. 4, right side). In both cases, the 
expected value of Di is 0, but instead Di mirrors D2, with positive values of Di that make it 
appear as though P2<^>Pi has also occurred. We can conclude both from theory and simulations 
that Di and D2 are inversely related, and therefore are non-specific in determining which pair of 
taxa have introgressed. 

The patterns counted by D12 could also be arrived at through introgressions other than the 
ones intended to be tested. For example, ABBBA is assumed by Eaton and Ree (2013) to be the 
result of a P?,=>P2 or Pa=>P2 transfer of the B state from a pre-introgression pattern of AABBA. 
However, ABBBA can also be formed from the patterns ABABA and ABBAA by Pi^>Pi> or P2=>P4, 
respectively. This means that for any given inter-group pair of taxa, introgressions in both 
directions will both raise or both lower the value of D12. This occurs because Eaton and Ree 
(2013) assume that the substitution A— >B occurred on the ancestral P34 branch, which also forces 
P3 and Pa to both be state B. This effectively collapses P3 and Pa to a single branch and makes 
D12 a four-taxon D-statistic. 

Simulations also confirmed these conclusions, showing that D12 values for a given pair of 
taxa are either positive or negative regardless of direction of introgression. For example, among 
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the four introgressions shown in Figure 4, P?,=>P\ and P$=>P\ exhibit negative values of Dn, 
respectively, as expected. However, P2=>P?, and Pi=>P\ are expected to be Dn ~ 0, but while they 
appear closer to zero, both cases still show a non-zero value of Dn (approximately at the same 
level as D\) with the same sign as introgression in the opposite direction. These simulations 
confirm that Dn will have the same sign for an introgression between given pairs of taxa 
regardless of direction. 

Discussion 

Robustness of the Dfoil Tests to Special Cases 

There are several special cases in which introgression occurs at particular points in the 
phylogeny and for which detection of these events might be more problematic. However, the 
triangulation approach of the Dfoil tests allows them to still function even in these special cases. 
If speciation of the two subgroups occurred at approximately the same time (i.e., Ti ~ T3), such 
that the difference between the two is not clearly discernible, the signs of all four statistics would 
remain unaffected. This would, however, negate the possibility of ancestral introgression, since 
both speciation events occur at approximately the same time. 

If introgression occurs very close to the second speciation event (i.e., r ~ T3) this would 
cause Dfi and Dot to reduce to 0 in the cases where introgression is P\IP2=$P$IPa. In this 
instance, Dfi and Dol are only reflecting the reality that the introgression involves a population 
of Pi or Pi that is so recently diverged that they are practically indistinguishable genetically. So 
as x approaches T3, the Dfoil statistics effectively convert to the signature of an ancestral 
introgression (+/+/ 0 / 0 or -/-/ 0/ 0). 
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One assumption that may be violated in real datasets is the expectation of uniform 
substitution rates across the tree. If a particular taxon has a much overall higher substitution rate 
than the others, then all distances involving that taxon would be relatively higher due to an 
increased rate of substitution rather than an introgression involving its sister taxon. There are 
two straightforward solutions to these lineage- specific effects. The first would be to weight the 
expected D-statistic values away from 50:50 when there is a prior expectation that substitution 
rates are not equal between sister taxa in a subgroup. The second would be to simply exclude the 
terminal-branch- substitution site patterns (AAABA, AABAA, ABAAA, and BAAAA) from all 
calculations. As long as these counts are excluded from both sides of the equation, the 
expectation of equality of the left and right terms of each Dfoil statistic is not violated. 

Ancestral Intro gressions 

Ancestral introgressions can also be detected by the Dfoil framework, though the 
direction of this introgression cannot be detected strictly by the Dfoil signature. When 
calculating a single, genome- wide estimate for each of the Dfoil tests, determining the direction 
of ancestral introgressions is not possible. However, if, for example, both Pn^Ps and P^Pn 
are present at different locations in the genome, the difference in their expected values for Dfo 
and Dil could be used to distinguish between these two cases in the different regions. 

More Taxa, More Models 

In order to detect introgression in a phylogeny of six or more taxa, the simplest option is 
simply to subsample a part of the tree in the appropriate configuration and either to use the four- 
taxon D- statistic or the five-taxon Dfoil tests. The Dfoil tests offer the added feature of 
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information about the direction of introgression in the subgroup of interest, when a five-taxon 
subset is available in the symmetric configuration. Aside from the subsampling option, we 
envision that a formal expansion of the model for more than five taxa would also be possible for 
certain tree topologies using the previously described principles of the D-statistic. The 
symmetric five-taxon tree represents a special case with respect to ILS, since the topology 
dictates that all discordant coalescences must occur in the root. This leads to a simple 
distribution of gene trees and provides a topology where each of the four in-group taxa can be 
compared in a straightforward manner with the two taxa in the opposite sub-pairs. Beyond this 
special case though, the probability distribution of gene trees becomes far more complex due to 
increasing topological constraints (Rosenberg 2002; Degnan and Salter 2005; Degnan and 
Rosenberg 2006). Therefore, while an explicit solution in larger phylogenies is not impossible, 
the more straightforward solution is simply to subsample four taxa of interest and an appropriate 
outgroup. 

Some care must be taken when subsampling from a larger tree, or when planning which 
taxa to sample from nature. Introgression from taxa not included in the sample — what Durand et 
al. (201 1) refer to as "ghost taxa" — could cause incorrect results using Dfoilox the D-statistic. 
When a pair of taxa introgress within the sampled clade, we infer introgression from closer 
relationships between these two taxa. However, if we only sampled the recipient of introgressed 
alleles, this taxon will appear to be unusually divergent from its sister taxon. Alternatively, a 
false positive of introgression could be detected if the ghost taxon has shared ancestry with one 
of the sampled taxa, or two taxa in the sample were recipients of introgressed alleles from the 
same ghost taxon. In general, however, we would expect that introgressions from taxa not in the 
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sample would simply result in a noisier signal of introgression or increased rate of false 
negatives. 

In the Dfoil tests, only biallelic site patterns are used to determine the phylogenetic 
relationships between taxa. This simple distance measure is particularly ideal for closely related 
species with few sequence differences and many biallelic sites. Since the site patterns used in 
Dfoil fundamentally derive from relative phylogenetic relationships and gene trees, these same 
underlying phylogenetic relationships could be used to build a Dfoil method using other models 
of sequence evolution. The conceptual framework would remain the same, only the model of 
inferring phylogenetic relationships between taxa would be altered. 

Conclusions 

In clades of closely related species, where ILS is prevalent or there are few sequence 
differences (or both), it can be difficult to detect introgression. The Dfoil system offers a simple 
means to infer introgression in a symmetric five-taxon clade, requires little computational power, 
and functions even with relatively little sequence divergence and high levels of ILS. When 
computed on multiple alignments of whole-chromosomes, the added spatial context will allow 
for the detection of localized introgressions throughout the genome. Spatial context also makes 
possible the detection of introgressions among different combinations of taxa at various locations 
throughout the genome. The Dfoil statistic is also designed specifically to detect the direction of 
introgression, when sequence data for sufficient taxa are available. 

The Dfoil tests can be used to test for general introgression across the genome in datasets 
without reference genomes, using RNA-Seq, RAD-Seq, or other targeted sequencing 
technologies. Since such loci may not have a known order along chromosomes, these data are 
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not suitable for locus-by-locus testing of introgression. For these data, a single genomic mean 
value for each of the four Dfoil tests can be computed (as was done by Green et al. 2010), and 
the average direction of introgression can also be determined. This implementation of Dfoil 
offers a more diffuse — but still informative — look at introgression from a smaller subset of 
data without spatial context. 

Supplementary Material 

Simulation scripts, simulated site pattern counts and D-statistics are available through 
Dryad (URL PENDING). The DFOIL script is available on BitBucket 
(http ://w ww .bitbucket.org/dfoil) 
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Table 1. The expected sign (+, — or 0) of Dfo, Dil, Dfi, and Dol for no introgression, eight 
inter-group introgressions, four ancestral introgressions, and two pairs of intra-group 
introgressions. Each inter-group introgression can be determined by a unique Dfoil "signature. 
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Figure Captions 

Figure 1. (a) A four-taxon phylogeny with three in-group taxa (P1-P3) and an out-group (O). (b) 
A five-taxon phylogeny with four in-group taxa (Pi— Pa) and an out- group (O). Ancestral 
branches are P12, P34, Pm, and P1234. Introgressions (dotted region) are shown from P^Pi in 
both (a) and (b). The time -before-present of the two or three speciations, respectively, are Ti, Ti, 
and T3, and time-since-introgression is z. Five-taxon trees can have (c) symmetric or (d) 
asymmetric topologies, shown with examples of allelic state patterns ABBAA and ABABA, 
respectively. 

Figure 2. (a) The three possible gene trees for a four-taxon phylogeny, including the gene tree 
concordant with the species phylogeny (i) and the two discordant gene trees (ii and iii). For each 
gene tree the biallelic site patterns used in the D-statistic are shown. Note that the "left" (+D) 
and "right" (— D) terms of the D-statistic sample gene trees whose relative probability is equal, 
(b) The fifteen possible gene trees for a five-taxon phylogeny, including the concordant gene tree 
(i) and discordant gene trees (ii-xv). The effects of discordant topologies for the Dfo statistic is 
shown as an example. Dfo samples sets of gene trees for the "left" (+Dfo) and "right" (-Dfo) 
terms, whose relative probabilities are equal regardless of their absolute probabilities. The 
relative values of gene tree pairs where ILS only occurs on a single branch (ii/iii and iv/v) 
depend on the relative length of the ancestral branches P12 and P34. 

Figure 3. The Dfoil statistics can detect both the taxa involved in and direction of introgression. 
(a) In a four-taxon symmetric topology, three classes of introgression can be defined: "intra- 
group" between taxa in the same subgroup, "inter-group" between taxa in different subgroups, 
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and "ancestral" involving the ancestral population of one subgroup, (b) Topology of tree used 
for simulations with species divergence times given in terms of millions of years ago (Ma) listed 
next to their respective nodes and inter-group introgression occurring between 0.20 and 0.25 Ma 
(dashed line), (c) The same topology showing the time of ancestral introgression between 0.53 
and 0.58 Ma. (d) Simulated chromosomes with n=500 10 kb regions for each of the sixteen 
possible introgressions. Between each introgression are n=100 simulated 10 kb sequences 
without introgression. Mean values over 20 consecutive regions for each Dfoil statistic are 
shown by solid and dashed lines. Each of the eight inter-group introgressions (first eight non- 
shaded columns from left) all show a unique Dfoil signature (i.e., a unique combination of +/-/0 
signs for Dfo, Dil, Dfi, and Dol) consistent with Table 1. The ancestral introgressions 
(involving Pu) can be distinguished between Pi2<=>P3 and Pi2<=>P4, but the direction cannot be 
determined by Dfoil signature alone. For the intra-group introgressions (rightmost four non- 
shaded columns) and no-introgression treatments (narrow, shaded regions throughout), all three 
Dfoil statistics are approximately 0 as expected. 

Figure 4. Partitioned D-statistics calculated for the same simulated chromosome and phylogeny 
as shown in Figure 3b. The Partitioned D-statistics suffer from a "mirror effect" due to the 
inverse relationship of Di and Di. For example, in the region under introgression Pi=>P3 
(leftmost column), Di is negative as expected, but Diis positive (and should =0). Also, Dn has 
the same sign for introgressions in both directions for each pair of taxa. 
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